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Abstract 

In the context of the phenomenon of Stochastic Resonance (SR) we study the correlation function, 
the signal-to-noise ratio (SNR) and the ratio of output over input SNR, i.e. the gain, which is 
associated to the nonlinear response of a bistable system driven by time-periodic forces and white 
Gaussian noise. These quantifiers for SR are evaluated using the techniques of Linear Response 
Theory (LRT) beyond the usually employed two-mode approximation scheme. We analytically 
demonstrate within such an extended LRT description that the gain can indeed not exceed unity. 
We implement an efficient algorithm, based on work by Greenside and Helfand (detailed in the 
Appendix), to integrate the driven Langevin equation over a wide range of parameter values. The 
predictions of LRT are carefully tested against the results obtained from numerical solutions of the 
corresponding Langevin equation over a wide range of parameter values. We further present an 
accurate procedure to evaluate the distinct contributions of the coherent and incoherent parts of the 
correlation function to the SNR and the gain. As a main result we show for subthreshold driving 
that both, the correlation function and the SNR can deviate substantially from the predictions of 
LRT and yet, the gain can be either larger or smaller than unity. In particular, we find that the 
gain can exceed unity in the strongly nonlinear regime which is characterized by weak noise and 
very slow multifrequency subthreshold input signals with a small duty cycle. This latter result is 



in agreement with recent analogue simulation results by Gingl et al. in Refs. [18, IS]. 
PACS numbers: 05.40.-a, 05.10.Gg, 02.50.-r 



I. INTRODUCTION 



Over the last twenty years or so, a large amount of work has been devoted to the study 
of the dynamics of noisy nonlinear systems driven by external periodic forces. One of the 
main reasons for this interest is related to the phenomenon of Stochastic Resonance (SR) 
n, |2|, |3|, |], H namely, the possibility of using the concerted action of noise and nonlinearity 
to augment selectively, for some parameter values, the output of the nonlinear system with 
respect to what it would be for a linear system dynamics. 

The two common quantifiers for Stochastic Resonance are the spectral amplification mea- 
sure P, H [7] and the signal-to-noise ratio (SNR) ^|, |8|. They are defined in terms of the 
Fourier components of the correlation function associated to the stochastic variable, x(t). 
Due to the periodicity of the driving force, the stochastic process x(t) is explicitly nonsta- 
tionary. Thus, the two-time function (x(t + T)x(t)) depends on both t and r. For very large 
values of t, this quantity is periodic in t with the period of the external driving. Thus, its 
cycle average over one period of t yields a function of just r: the correlation function, C(r). 
The analysis of its structure reveals that C(r) is the sum of two terms 0: One term is 
periodic in r with the same period as the driving force and it is called the coherent part, 
C C oh(j). The other term, the incoherent part Ci nco h(r), decays to zero for r — > oo. The SNR 
of the output process, x(t), is defined as the ratio of the amplitude of the Fourier mode of the 
coherent part at the driving frequency, and the power spectral density of the incoherent part 
taken also at the driving frequency. By definition, the SNR is thus a dimensional quantity. 

The SNR of an input signal (SNR| inp ) containing the sum of the external driving and the 
Gaussian white noise can easily be evaluated. A convenient dimensionless parameter can 
then be defined as follows: the gain, G, defined by the ratio of SNR of the output over the 
(SNR|j np ). For the case that the Langevin dynamics is linear in x driven by additive white 
Gaussian noise, then the output SNR is exactly the same as SNR|j np ; i.e. the gain assumes 
precisely the value unity. In a general nonlinear case, neither the SNR nor the gain can 
be evaluated exactly by analytical means. As a consequence, their evaluation necessarily 
requires approximate procedures. 

It was pointed out previously that the gain of a noisy nonlinear dynamical system subject 



to subthreshold signals can not exceed 1 R 10 1. This feature has been rationalized using the 



ideas of Linear Response Theory (LRT), thought to be valid for weak driving amplitudes and 
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not too small noise strengths. It should be pointed out, however, that the validity of LRT 
critically depends also on the value of the frequency, as has convincingly been demonstrated 



in recent works Rill 12 



In the context of LRT theory it has been pointed out in || that a corollary of LRT is 
that "for small amplitude signals, the signal-to-noise ratio at the output of a system driven 
by a stationary Gaussian noise does not exceed that at the input, even if the system displays 



SR". Moreover, in Ref. jKj, the authors state that for "small signal in a Gaussian noise 
background, it is a theorem that the SNR at the output of a non-linear device must be less 
than or equal to the SNR at the input" . On the other hand, studies on nondynamical systems 



fl3| , |14| , jig , jig] , on dynamical systems driven by large amplitude sinusoidal forces ||17| , and 
on dynamical systems driven by pulsed (multifrequency) periodic forces with subthreshold 



amplitudes [18. 19|, have reported gains larger than unity. Clearly, for this to occur, the 
stochastic system must operate in a regime where LRT does not apply. It is therefore 
of interest to delineate carefully the limit of applicability of the LRT description of the 
correlation function, the SNR and the gain of a nonlinear noisy driven system. 

In this paper we have tackled this challenge by carrying out a detailed numerical eval- 
uation of the correlation function C(t) and its components, C co h{r) and Cj nco / l (r), of the 
SNR and the gain of a bistable noisy system which is driven by time periodic forces. The 
numerical predictions have been compared with those provided by the LRT approximation 
which accounts for the full spectrum of all relaxation modes. 

As it is well known, LRT requires the knowledge of the system susceptibility, or al- 
ternatively, of the correlation function of the noisy system in the absence of driving, 
K(t) |2|, [11], 0, ^T], |22|. None of these quantities are known exactly for nonlinear 
systems. For sufficiently small values of the noise strength, suitable analytical approxima- 
tions to K(t) can been used [|[ [11], [L^, |20], [21| . On the other hand, for large values of the 
noise intensity, we have evaluated K(t) from the numerical solution of the Fokker-Planck 



equation using an adaptation of the split operator technique of Feit et al. p3fl , as it has 
been detailed in [p4}| . In this paper, we also present a detailed proof of the statement that 
within LRT, the gain Qlrt — 1 ; by use of the full spectral approach; this proof differs 
from alternative attempts in Refs. || [K| which use additional restrictions such as a linear 
response theory for the fluctuations themselves. 

The "typical" procedure to evaluate the SNR involves the Fourier analysis of a very long 



3 



record of the stochastic trajectory, x(t). Using the Fast Fourier Transform (FFT) of the 
record, the corresponding periodogram is constructed. There are several drawbacks with 
this procedure. There are subtleties inherent to the interpretation and evaluation of the 
periodogram (see for instance the critical comments in [2S|). There are also major problems 
associated with the fact that the power spectrum contains 5-peaks at the driving frequency 
and its higher harmonics arising from the coherent part of the correlation function. The 
contribution of the incoherent part at those frequencies is embedded in those peaks, and 
it is not a simple task to estimate the separate contribution to the peaks of the coherent 
and incoherent parts of the periodogram. The evaluation of the SNR-gain requires a good 
knowledge of both contributions, and any small error in the estimation of the incoherent 
contribution, yields unreasonable values for the gain. Indeed, in our opinion, a much better 
estimate would be obtained if the periodic part of the output signal were subtracted from 
the data before performing its FFT. 

In this work we propose such an alternative procedure. The Langevin equation is nu- 
merically integrated for a large number of noise realizations. The time evolution of the 
correlation function and its coherent part are directly evaluated from the numerical solution 
after averaging over the noise realizations. The incoherent part is obtained from the differ- 
ence Cincohij) = C{r) — C co h{j). As the definition of SNR requires just the amplitude of 
the Fourier mode of C co h{r) and the spectral density of C inco h{r) at the driving frequency, 
the SNR can readily be evaluated with two numerical quadratures; i.e. there is no need to 
construct the full spectrum. 

The paper is organized as follows. In the next Section we introduce the model and 
provide definitions of the quantities of interest. In Section III, the main points of the LRT 
description of the correlation functions are detailed. We also present in this Section a novel 
and straightforward proof of the fact that Qlrt < 1> based on the spectral properties of 
the Fokker-Planck operator, and its adjoint, in the absence of driving. In Section IV, we 
present the numerical procedure used to obtain the correlation function, the SNR and the 
gain from the numerical solution of the Langevin equation. The very efficient algorithm 
used in this work is summarized in the Appendix. The numerical results are compared with 
the predictions of LRT for a variety of parameters and two distinct types of driving forces: 
a monochromatic force and a periodic sequence of pulses. Finally, we present conclusions 
for the main findings of our work. 
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II. CORRELATION FUNCTION, SIGNAL-TO-NOISE RATIO AND GAIN 



Let us consider a system characterized by a single degree of freedom, x, subject to the 
action of a zero average Gaussian white noise with (£(£)£(s)) = 2D5(t — s) and driven by 
an external periodic signal F(t) with period T. In the Langevin description, its dynamics is 
generated by the equation 

x(t) = -U'[x{t)]+F(t)+£{t). (1) 

The corresponding linear Fokker-Planck equation (FPE) for the probability density P(x,t) 
reads 

^P(x,t) = l{t)P{x,t), (2) 



where 

p> r a i 

(3) 



In the expressions above, U'(x) represents the derivative of the potential U(x). The period- 
icity of the external driving F(t) allows its Fourier series expansion in the harmonics of the 
fundamental frequency Q = 2n/T, i.e., 

oo 

F (t) = E Ifn cos(nnt) + g n sin(nfit)] , (4) 

n=l 

with the Fourier coefficients, f n and g n , given by 

f n = % [ T dtF(t)cos(nnt) 
1 Jo 

9n = T L dtF ( t ) sin ( nnt )- ( 5 ) 

Here, we are assuming that the cycle average of the external driving over its period equals 
zero. 

The two-time correlation function (x(t + r)x(t))oo in the limit t — > oo is given by 

/oo /*oo 
dx'x'P oc (x', t) I dxxPx\x(x, t + t\x', t), (6) 
-oo J — oo 

where Pqo(x, t) is the time-periodic, asymptotic long time solution of the FPE and the 
quantity P 1 | 1 (x,t + r\x',t) denotes the two-time conditional probability density that the 
stochastic variable will have a value near x at time t + r if its value at time t was exactly 
x' . It can been shown [|2|, |7H that, in the limit t — > oo, the two-time correlation function 



(x(t + r)x(t))oo becomes a periodic function of t with the period of the external driving. 
Then, we define the one-time correlation function, C(r), as the average of the two-time 
correlation function over a period of the external driving, i.e., 

C(r) = ^J T dt(x(t + r)x(t)) oo . (7) 

The correlation function C(r) can be written exactly as the sum of two contributions: a 
coherent part, C co h{T), which is periodic in r with period T, and an incoherent part which 
decays to for large r. The coherent part C co h{r) is given by p], [7| 

C coh {r) = l- / T dt(x(t + r)) 00 (a;(t)) 00 , (8) 
1 Jo 

where (x(t))oo is the average value evaluated with the asymptotic form of the probability 
density Poo(x, t)- 

It is possible to carry out a formal analysis of C(r) and its coherent and incoherent com- 
ponents by making use of the spectral analysis of the Floquet operator associated with the 
Fokker-Planck dynamics. But an explicit evaluation of the correlation function is generally 
impossible; thus, one has to rely on numerical results obtained from integrating either the 
Langevin or the FPE, or by use of approximate analytical descriptions. 

According to McNamara and Wiesenfeld [H] , the SNR is defined in terms of the Fourier 
transform of the coherent and incoherent parts of C(r). As the correlation function is even 
in time and we evaluate its time dependence for r > 0, it is convenient to use its Fourier 
cosine transform, defined as 

C(uo) = - drC(r)cos(cjr); C(r) = / dcu C(u) cos(cjt). (9) 

7T JO JO 

For the SNR we define: 

SNR = ^o+^d^M (1Q) 

Note that this definition of the SNR differs by a factor 2, stemming from the same con- 
tribution at u = —Q, from the definitions used in earlier works |2|, 0. The periodicity of 
the coherent part gives rise to delta peaks in the spectrum. Thus, the only contribution to 
the numerator in Eq. (|10D stems from the coherent part of the correlation function. The 
evaluation of the SNR requires the knowledge of the Fourier components of C co h{r) and 
Ci n coh{ T ) at the fundamental frequency of the driving force. The entire Fourier spectrum 
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is not needed. The evaluation of the SNR thus requires the evaluation of two well defined 
numerical quadratures only; i.e., 

c N d _ I So dT C coh(r) cos(^r) 

^ J °° dr C incoh (r) cos(nrY 1 ) 

The signal-to-noise ratio for an input signal F(t) + £(£), is given by 

SNR\ mp = <f \ + / l) . (12) 

The so-called gain is defined as the ratio of the SNR of the output over the SNR of the 
input; namely, 

SNR\ 

inp 

III. LINEAR RESPONSE THEORY BEYOND THE TWO-MODE APPROXIMA- 
TION 

The Linear Response Theory provides a general procedure to describe the correlation 
function in an approximate way. The basic quantity of LRT is the system response function, 
x(t)- It is related to the equilibrium time correlation function of the system in the absence 
of external driving, K(t), via the fluctuation-dissipation theorem (FDT) || |7|, |20], |21| 



i.e. 



*(«> = < ° : i£ ° (") 

-±K(t) : t > 
The equilibrium time correlation function K(t) is defined as 

dx'x'P (eq) (x') / dxxP$(x,t\x'), (15) 

where P^ eq \x) is the equilibrium distribution of the non-driven system, 

P^\x) = Ne-^, (16) 

and P$(x,t\x') is the conditional probability density to find, in the absence of driving, the 
variable near x at time t, if it was initially at exactly x'. Here we are assuming that the 
potential U(x) is even in x, so that (x) eq = 0. 

Within LRT, the long time average value (x(t))^ RT ^ is given by 

/•oo 

(x(t))^= / dr X (r)F(t-r). (17) 
Jo 
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Insertion of the Fourier expansion Eq. (|4]) into Eq. (|17|) leads to 

oo 

«t))£ RT) = E [M { n RT) cos(nfit) + Ni LRT ^ sin(nfit) 

n=l 

where the coefficients MfJ jRjn> and N^ LRT ' are given by 

M (LRT) = f (r) _ (i). jv(£«T) = f (i) , (r) 
1H n JnAn ilriAn i iv n J nXn ~ tin An • 



(18) 



(19) 



In these formulas, we have introduced the quantities x« "* an d Xn"* defined as 

Xi r) = /°°drx(r)cos(nfir) (20) 



dr x( r ) sin(nfh~). 



(21) 



The use of the FDT in the above expressions allows us to write immediately 

(r) (x 2 ) eq - nQ / °° dt K(t) smjnttt) 



Xn 
An 



fl^l poo 

— / dtK(t) cos( nftt). 
D Jo 



(22) 
(23) 



It then follows from Eq. (H) that within LRT, the coherent part of the correlation function 
is given by, 



n=l 



cos(nfh~). 



(24) 



As discussed in Refs. [§, 0, |2(| , LRT amounts to keeping the leading term in the perturbation 
treatment of the dynamics of the stochastic process x(t) in powers of the driving amplitude. 
Then, within the spirit of perturbation theory, the leading term in the expansion of the 
incoherent part corresponds to the correlation function of the system in the absence of 
driving force, i.e., C\^(t) = K(t). 

Taking into account that C^ RT (r) is periodic in r, it follows from Eqs. (|) and (p4[) that 

> ( U 1 ' ' '""•' ,-> 



1 OO r 

n=l L 



[6(nQ - u) + 5(nfi + cj)] . 



(25) 



Thus, it follows from the definition of the SNR, Eq. (|T0D , that, within LRT, we have 

\m[ lrt) ) 2 + (N[ LRT) ^ 



Xi 



(r) 



X 



(0 



AD X f 



(26) 



where K(Q) is the Fourier cosine transform of K(t), defined according to Eq. (^|). In arriving 
at Eq. @ we have also used Eqs. flT^, (0) and (0). 

Taking into account Eqs. (P3)> (0) an d (|26|), one readily finds that the gain within LRT 
is given by 



q(LRT) 



SNR (lrt) n 



Xi 



(r) 



Xi 



(0 



Xi 



(27) 



SNR\ inp 

This is a general expression for G^ Li?T ^ valid for any shape of periodic driving. 

The last expression will allow us to show that G^ iRT ^ can, indeed, not exceed unity. 
Although this assertion has been discussed previously in Ref. || [10| , we next will present a 
detailed and hopefully very clear proof for this prominent assertion. 



As shown in the Appendix of Ref. 0, see also in Ref. JTT|, |I2| , the susceptibility, x(*)> 
can be expressed as 

oo o 

xW = -E e " v (°l^b)(pl|-l )' ( 28 ) 

p=i ax 



where \p) = ip p (x), (p\ = ipUx) and \ are the eigenfunctions and eigenvalues of the FP 
operator Cq associated to the undriven dynamics and its adjoint, £q, i.e., 



C ip p {x) = -\ p ip p (x), CoipUx) = -Xpiphx). 



(29) 



Using the above representation of the susceptibility in Eqs. (|20|) and (|2~T|) with n — 1, we 
find 



(r) 
Xl = 


oo 

-E 




d 

{0\x\p){p\ — 




(i) 
Xl = 


oo 

-E 

p=l 




d 

(0\x\p){p\ — 


\ 2 p + W 



A, 



(0|x|p>(p|^|0> 
(0|xb)(p|^|0) 



(30) 
(31) 



Here, we have used the inequality 



Np)(p|^|0) <o, 



(32) 



which can be proved as follows. Multiplying the first equation in Eq. (|29]) by x and carrying 
out an integration by parts, one obtains 



X p (0\x\p) 



dx x CqiP p (x) 
d 



dx U' (x)i/jo(x)i/jI(x) 



Jx^l{x)-Ux)=D(p\-\Q) 



(33) 
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where we have taken into account that ipo(x) = P^ eq \x) and ip p (x) = ■ip (x)ipj ) (x) , so that 
tpl(x) = 1. Therefore, (0\x\p) (p\£\0) = -X p ((0\x\p)) 2 /D < 0. Using in Eqs. © and (P) 
the Cauchy-Schwarz inequality, we find 



(xf 



r) 



00 A„ 

E- 

P =i 



<0|*|p>(p|£|0> 



1/2 



A2 + fi2 



< g a;|<o|x|p)<p|^|o)| g 



P =i 



E 

P =i 



A2 + fi2 



n|(o|x|p>(p|^|o> 



9=1 
1/2 



A2 + fi2 



< E 



Q 2 ko|s|p><p|&|0> 



A 2 + n 2 



E 

9=1 



(0|x|p)(p|^|0) 



M?X?I^I°> 



(0|x|p)(p|^|0) 



(okl?X?l^|o> 



1/2 



(34) 



1/2 



(35) 



Taking into account that (0|x|0) 



the completeness relation yields 



E*i (0\x\q)(q\£\0) 



E~ (0k|g)(g||:|0) = -(0\x£\0) = 1. Thus, by adding 



Eq. to Eq. (g5|), one obtains 



, {r)V ( m 2 - \{0\x\p)(p\f x \0) 

bci) +{xi) < E 1 
P =i 



(i) 

Xi 

A2 + fi 2 Q 



(36) 



Finally, inserting Eq. (|36|) into ([27]), we obtain the seminal inequality that G^ LRT ^ < 1. 

Put differently, the gain of a nonlinear system operating in a regime where LRT provides 
a valid description cannot reach values greater than 1. This result is valid for any periodic 
external driving. Notice that this finding does not preclude the possibility of obtaining 
values for the SNR-gain larger than unity when the conditions are such that the use of LRT 
is not sensible. 

IV. NUMERICAL RESULTS 



In this Section, we will carry out the numerical evaluation of the different magnitudes 
defined above. Our goal is to compare the predictions of LRT with the results obtained from 
the numerical solution of the Langevin equation, Eq. ([!]). We will consider the dynamics in 
the bistable potential U(x) = —x 2 /2 + x 4 /4 driven by time periodic forces. 

The evaluation of the different magnitudes using LRT requires the knowledge of K(t) [cf. 
Eqs. ( |19| , p4] , ^6 , $7\)]. For nonlinear problems, explicit expressions for K(t) are unknown, 
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but useful approximations have been presented in the literature. For the bistable potential, 
U(x) = —x 2 /2 + a; 4 /4, Jung and Hanggi [[TL]] have used the two-mode approximation. It is 
based on the existence of a large difference in the time scales associated to inter-well and 
intra- well motions, and it is expected to be valid for small values of the noise strength D. 
With this model, one finds 

K(t) = g x e- x ^ + g 2 e- aT (37) 



where 



A lRj ^(l-|D)exp(-l/(4D)), (38) 

7T Z 



and a = 2. The weights, g\ and g 2 , can be obtained from the moments of the equilibrium 
distribution in the absence of driving using the expressions 

_ ) eq (x ) eq — (x ) eq (qc\\ 

92 — -\ 1 : [dy) 

Ai — a Ai — a 

9i = (x 2 ) eq - 92- (40) 

To leading order in D, we can replace Ai by A^ = V^/n exp(— 1/ (4D)), g\ « 1 and g 2 ~ -D/«. 
This is the limit considered in Ref. PB[ . In the results reported below, we have also 
considered values of D so large that the two-mode approximation becomes inadequate. 
Therefore, the correlation function in the absence of driving has been evaluated numerically 
from the FPE in the absence of driving following the procedure discussed in Ref. |^4| . 



The numerical evaluation of the correlation function C(r) and its coherent and incoherent 
parts proceeds as follows. Stochastic trajectories, x^'(t), are generated by numerically 
integrating the Langevin equation for every realization j of the white noise £(£), starting 
from a given initial condition xq. The numerical solution is based on the algorithm developed 



by Greenside and Helfand [p7]> p8| . The essence of the algorithm is briefly sketched in the 
Appendix. After allowing for a relaxation transient stage, we start recording the time 
evolution of each random trajectory for many different trajectories. Then, we construct the 
two-time (t and r) correlation function, i.e., 



1 N 

(s(t + T)s(t))o„ = --£xM(t + T)xW(t), (41) 
iV 3=1 



as well as the product of the averages 



1 N \ I 1 N 



(x(t + T)) 00 (x(t)) 00 =[i-Y,x (j) (t + T)\ hr£* W (*)l, ( 42 ) 



3=1 J \ 3=1 
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where N is the number of stochastic trajectories considered. The correlation function C{t) 
and its coherent part C co / l (r) are then obtained using their definitions in Eqs. (0) and 
performing the cycle average over one period of t. The difference between the values of 
C(t) and C co h(r) allows us to obtain the values for Ci nco h(r). It is then straightforward 
to evaluate the Fourier component of C co /i(r) and the Fourier transform of Ci nco h(T) at the 
driving frequency by numerical quadrature. With that information, the numerator and the 
denominator in the expression (|ToD for the SNR and the gain (Eq. fll3"D ) are obtained. 

We shall analyze two different types of periodic driving forces. First, let us consider 
the well known situation with a monochromatic, single frequency force, Acos(Ot), with 
amplitude strength A and angular frequency Q @. In this case, the formulas in Sect. Ill 
simplify considerably because fi = A, while all the other Fourier components of the driving 
force vanish. The second case corresponds to a periodic force with period T, with a sequence 
of pulses of length t c < T/2 ; namely, 



F(t) 



A : < t < t c 

-A : Z< t< T +tc (43) 
: otherwise. 



In this case, we have 



2A 2A 
j\ = — sm(Qt c ); gi = — [1 - cos(^t c )] , (44) 

7T 7T 

where Q = 2tt/T is the fundamental frequency. This force is characterized by its amplitude, 
its period and its duty cycle, which is defined as 2t c /T. Recently, Gingl et al. [|TH], [T^ 



have carried out analogue simulations of systems that are subjected to wideband Gaussian 
noise and driving forces of this second type. They report values for the gain which greatly 
exceeds unity, for driving amplitudes below its threshold value. If this is the case, then 
strong deviations from the LRT should be observed as well. 



A. Monochromatic driving 

In Fig. |I] we depict the results obtained for a monochromatic driving force with angular 
frequency Q = 0.1, noise strength D = 0.2 and several values of the amplitude. In the 
deterministic dynamics (D = 0), an external periodic force with the indicated frequency 
induces sustained oscillations between the minima of the potential for A > A t h ~ 0.419. 

12 



Note that this nonadiabatic frequency raises the threshold value for superthreshold driving 
beyond its adiabatic lower limit of A^f 1 = ^4/27 ~ 0.3849. Thus, we will take this value 
as the amplitude threshold value at the frequency Q — 0.1. In panel (a), we plot the 
numerators of the SNR in Eq. ( JlQ| ) and of SNR( LiJT ) in Eq. fl26|) vs. A 2 . The solid straight 
line represents the LRT result, while the circles correspond to the numerical results. The 
graph reveals that for amplitude strengths A < 0.1 the predictions of LRT match well the 
numerical results, as can be expected. When the amplitude increases, the deviations of LRT 
from the precise numerical results are large. LRT predicts a much larger amplification of the 
output amplitude than the one obtained numerically. In panel (b), we plot the denominators 
of the SNR in Eq. (0) (circles) and of SNR( Li?T ) in Eq. © (solid line) vs. A 2 . In LRT, the 
denominator is independent of A. Once again, the predictions of LRT match the numerical 
results for A < 0.1. For larger values of A, the influence of the driving amplitude on the 
relaxation of Ci nco h(r) is very strong and the numerical results for the denominator are 
much smaller than the ones obtained within LRT. It is then clear that LRT will yield a 
valid description of the signal-to-noise ratio for small driving amplitudes only as depicted in 
panel (c). We notice that the values for the SNR provided by the numerics are larger than 
SNR( Li?T ) . This is so although linear response theory predicts larger spectral amplifications, 
see in |7[, of the average output than what really occurs. The modifications in the behavior 
of the incoherent part of the correlation function with respect to its behavior in the absence 
of driving are more than enough to compensate for the behavior of the numerators. In 
panel (d) we plot the gain vs. A 2 . There exists an optimum value for the driver amplitude 
(A ~ 0.8) at which the gain becomes maximized. Nonetheless, the gain is always smaller 
than unity. LRT requires that G^ LRT ^ < 1. These strong deviations of the predictions of 
LRT about the behavior of the two components of the correlation function with respect to 
the numerical results tell us that LRT cannot be invoked to explain the fact that the gain is 
smaller than 1 for the range of parameter values considered in this figure; i.e. a gain below 
1 occurs here within the nonlinear regime. 

In Figs. |2| and |3| we analyze the same quantities as in Fig. [1], but now for larger noise 
values, D = 0.6 and D = 1.0, respectively. The most important difference with respect to 
the plots in Fig. [1] is that for these larger values of the noise, the gain can exceed unity for 
values of the amplitude well above its threshold value. This superthreshold feature has been 
corroborated already in Ref. [|T7[]; a gain above 1 seemingly does not occur for monochromatic 
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subthreshold driving. 



B. Pulsed, multichromatic periodic driving 

Next, we proceed to consider the case of pulsed driving forces. In Figs. |] and f| we compare 
the dependence of the output on the driving amplitude as given by the LRT approximation 
with the numerical precise results. The system is forced by a multifrequency driver with a 
period T = 27r/0.1 ~ 63 and a duty cycle of 10%. As in the case of a single frequency driving, 
the values of the different quantities obtained from the numerics deviate significantly from 
the predictions of LRT as the amplitude of the driver is increased. Nevertheless, perhaps 
the most relevant difference with respect to the monochromatic case is that we again do not 
find gains larger than 1 in the range of parameter values considered in these figures. 



C. The case of strong nonlinearity 

A particularly interesting situation arises in the analogue studies of pulsed driving forces 



with a very small fundamental frequency: in Refs. [0, [19| Gingl et al. report gains that 
significantly exceed the value 1 for a subthreshold, multifrequency driving force of very large 
period T = 27r/0.0024 ~ 2618 and a small duty cycle of 10%. This large gain is accompanied 
by a non-monotonic behavior of the SNR with the noise strength D. Therefore, this situation 
must correspond to a very sensible discrepancy of the actual behavior with respect to the 
LRT predictions. We have carried out detailed and careful numerics of the Langevin equation 
in this extreme regime for such a driving force with a subthreshold amplitude A = 0.35 
and a noise strength D = 0.02. With the parameters considered, the problem becomes 
computationally very demanding indeed: this is so because of the very large period of the 
driving force. Moreover, in order to obtain reliable numerical results for the incoherent part 
of the correlation function a large number of stochastic trajectories needs to be generated. 
Our findings are summarized with the Table |I[ 

To obtain a reliable convergence of the corresponding SR quantifiers at least up to 50000 
random trajectories need to be considered. A smaller sampling size can induce severe errors, 
see in Table I. The main result is a numerically evaluated gain of 8.62; in clear contrast, 
the result predicted by LRT is the very small value of 0.018; i.e. LRT strikingly fails, cf. 
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TABLE I: 





trajectories 


numerator 


denominator 


SNR 


eain 

o 


numerics 


1000 


0.78 


0.33 


2.32 


12.16 




5000 


0.78 


0.35 


2.26 


11.84 




10000 


0.78 


0.47 


1.67 


8.77 




50000 


0.78 


0.48 


1.65 


8.62 


LRT 




0.00061 


0.177 


0.0034 


0.018 



in Table I for the corresponding values of SNR and its constituents. The SNR value of 



the analogue simulation in Refs. [18|, [19| carried out with a pulsed input signal with the 
same characteristics as the one considered here, and wideband Gaussian noise with a related 
strength roughly similar to ours, yields an experimentally determined gain of ca. 19, cf. 
Fig. 4 in fT9"l| . This value is again significantly larger than 1 and compares favorably with 
our results in Table I. Note, however, that the sampling size of ca. 1000 realizations used in 



Refs. [|TI| [L9| has been chosen substantially smaller than the number of realizations needed 



to achieve good numerical convergence, cf. Table |; this in turn may explain the overshoot 
of the experimentally determined gain value. 



V. CONCLUSIONS 



Let us summarize the main results of this work: (i) First, we have provided an analytical 
proof based on LRT beyond the commonly employed two-mode approximation that the gain 
of a noisy, periodically driven nonlinear system which operates within the regime of validity of 
LRT cannot exceed unity. This result holds for arbitrary noise strength D and is independent 
of the shape of the input signal, (ii) We have implemented a very efficient algorithm due 



to Greenside and Helfand [£7| [28[ to numerically integrate the Langevin equation. From 
the numerical solution, we have evaluated the time evolution of the correlation function 
and its coherent and incoherent components, (iii) We have also put forward a procedure, 
alternative to the usual one, to calculate the SNR. The numerator and denominator of the 
SNR are calculated by use of only two numerical quadratures, (iv) A detailed comparison 
between the predictions of LRT and the numerical results have been carried out. We have 
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assessed regions of parameter values where LRT gives an erroneous description, yet the gain, 
nevertheless, is less than unity. On the other hand, there exist regions in parameter space 
where the gain indeed exceeds 1 if driven with a superthreshold amplitude strength; this 
finding is in agreement with prior results in Ref. jr7| . These regions are again characterized 
by substantial deviations from LRT. 

Moreover, as previously established by use of analogue simulations by Gingl et al. |TB|, [H3 



we also find the surprising result, valid for dynamical systems, that SNR-gains larger than 
unity can indeed occur for subthreshold (!) polychromatic input signals: For this feature to 
occur one seemingly needs, however, weak noise and a slow periodic driving signal with a 
very small duty cycle. In this context, the necessity of a sufficiently large number of sampling 
trajectories in order to obtain reliable, convergent results has also been stressed. It is in this 
very regime of small frequency driving and weak noise where the LRT description indeed 
fails notably . 



APPENDIX A: THE METHOD OF GREENSIDE AND HELFAND 

The procedure proposed by Greenside and Helfand for numerically integrating stochastic 
differential equations has been discussed in detail by their authors in |27]. 59] . For the sake of 



completeness, we will briefly sketch in this Appendix the main reasoning of their procedure. 
By analogy with deterministic Runge-Kutta algorithms, Greenside and Helfand developed 
schemes to estimate the value of the stochastic variable at time t + h if its value at time t 
is known. This is achieved by evaluating the right hand side of the Langevin equation at 
selected points within each interval of length h, so that, all moments of x(t + h) — x(t) are 
correct to order h k . 

As our Langevin equation contains an explicit time dependent driving force, it is conve- 
nient to rewrite it as a two-dimensional problem with variables (2/1, 2/2) = V, where y\ — x 
and t/2 = t. The Langevin equation, Eq. fll]) is then written in vector form as 

f t =G(y) + m (Ai) 

where G = (G h G 2 ) = (-U' (x) + F(t), 1) and E(t) = (£(*), 0). 
The formal solution of Eq. ( |A1|) yields 

y K (h) = y K (0)+ [ h dsG K (y(s))+w^(h), (k=1,2) (A2) 
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with 



w<?{h)= f h dsE K (s) (A3) 
J o 

The right hand side of Eq. (|A2|) can be expanded as 



y K (h) = y K (0) + hG K (y(0)) + h 2 £ d 3^1 G ^ m) + . . . + s K {h) (A4) 

The last term, S K (h), represents the stochastic part. It is a series in h 1 / 2 with the order of 
the terms determined in probability. 

By analogy with the Runge-Kutta procedures for deterministic differential equations, 
Greenside and Helfand propose an /-stage algorithm to write the solution of Eq. (|Al|) as 



y*(h) = y*(0) + h(A l9lK + ... + A l9lK ) + h*Z%Y 0K (A5) 



with 



9u = G K (jj/ M (0) + h^lY^ 
92k = G K ({^(0) + hfingxp + h^Y 2 ^ 

9in = G R Nj/ M (0) + hPngip + ... + hfr^g^ + h^Z^Y^ (A6) 

Here, {{y^}) is the set (x,t). The Y^ are Gaussian stochastic variables with zero average 
which are numerically generated by writing 

m 

(A7) 

3=1 

where Zj K are m independent Gaussian random variables of zero average and unit variance. 
The parameters A i: and Ajj appearing in Eqs. ( P£5[ - |A7[ ) are independent of the component 
index k. They are obtained by expanding Eq. (|A5|) to the desired order h k . This expansion 
gives rise to a deterministic and a stochastic part, S K . Equating the coefficients of this 
expansion with those of the deterministic part in Eq. ( |AJ ) leads to a set of equations for 



the parameters A+, and A^-. Further equations are obtained by equating the moments of 
(S™) with those of the stochastic part in the expansion in Eq. (|A4[) 

A procedure correct to order h k in the step size h, involving /-stages and m Gaussian 
independent variables is termed a kols m G algorithm. In this paper we have integrated the 
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TABLE II: Parameter values given by Greenside and Helfand [28| for their Sq^s^g algorithm 



Ai 


u.u 


A 2 


r\ O A A A c 


A 

A3 


U.ly44oU 


A 

A4 


U.lolUoZ 


ft 

P21 


(\ ~ \ r"-7-\ n 


ft 

P31 


-U.oy7o(JU 


Q 

P32 


u.4z /oyu 


ft 

PA1 


-l.OO MOl 


ft 

P42 


1.417263 


ft 

PA3 


1.170459 


\ 

A01 


1 n 
l.U 


^02 


n n 
U.U 


An 


0.0 


Al2 


0.271608 


A21 


0.516719 


A22 


0.499720 


A31 


0.030390 


A32 


-0.171658 


A41 


1.0 


A42 


0.0 



Langevin equation using a ?>o^s^g algorithm with the values for A i} /3y and A^- given in 
Table |I| taken from |28]. With this choice of parameters, the deterministic part is of order 
/i 4 , as in the fourth order Runge-Kutta procedure for ordinary differential equations. 
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FIG. 1: The dependence of several SR quantifiers versus the square of the driving amplitude, A 2 , 
given by LRT (solid line) and by the numerical solution of the Langevin equation (circles) . In panels 
(a) and (b), we plot, respectively, the numerator and denominator appearing in the definition of 
SNR, cf. Eq. (|TT|). The behaviors of the SNR and the gain are depicted in panels (c) and (d). 
The driving force is monochromatic with frequency = 0.1 and the white noise strength is kept 
constant at the value D = 0.2. In all panels, the vertical dashed line indicates the square of the 
value of the dynamical threshold amplitude, Ath, at the angular driving frequency, £1. In panel 
(d), a dotted horizontal line is drawn at the gain value of 1 as a guide to the eye. 
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FIG. 4: The same as in Fig. |], for the case of a pulsed, driving force with period T ~ 63, duty 
cycle 2t c /T = 0.1 cf. Eq. (|43|), and a noise strength D = 0.6. 



C/) 



0.05 
0.04 
0.03 
0.02 
0.01 


0.06 
0.04 
0.02 



(a) 










o - 




so 












(c) 










/ o - 









LU 

a 



(b) 


W o 

o 










o 



(d) 



0.5 1 1.5 2 0.5 1 1.5 2 



0.8 

0.7 

0.6 
0.9 

0.8 

0.7 

0.6 



FIG. 5: The same as in Fig. [l] for a pulsed driving force with period T ~ 63, duty cycle 2t c /T = 0.1, 
cf. Eq. (HH), and a noise strength D = 1.0. 
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